!
!  AutoFlat     Mean of the scaled set of scientific images.
!  Copyright (C) 1997-x  Filip Hroch, Masaryk University, Brno, CZ
!
!  This file is part of Munipack.
!
!  Munipack is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation, either version 3 of the License, or
!  (at your option) any later version.
!  
!  Munipack is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!  
!  You should have received a copy of the GNU General Public License
!  along with Munipack.  If not, see <http://www.gnu.org/licenses/>.

Program AutoFlat

! Experimental version #5.
!
! Experimentalni verze 5.
! predchozi verze asi moc nefungovaly!!!
!
! Toto je nejdrsnejsi verze. Pocita primo median a tak
! potrebuje opravdu hodne pameti -- typicky desitky Mb.
!

! Vyvoreni autoflat-fieldu z mnoha CCD snimku.
!
! FH, Rijen 1997
! FH, duben 1998
! FH, cerven 1998
! FH, zari 1998
! FH, rijen 1998, 7ho - vypada to velmi dobre, i drsne testovaci snimky 
!                   to uspesne projelo
! FH, cerven 1999 - dynamicke alokovani 
! FH, prosinec 2004 - hlavicka
!

  use LibList
!  use optfile
  use input
  use stat
  use robustmean

  implicit none

  !    Version:
  Character(len=*),parameter :: version = 'mflat Version 0.3, (C)1997-2010 F.Hroch (hroch@physics.muni.cz)' 

  ! Initial parameters. 
  !    Output image name:
  Character(len=32768) :: nameflat='mflat.fits'

  ! No. of image dimensions
  integer, parameter :: DIM = 2

  !    Output mean:
  Real :: level = 1.0
  ! LEVEL je vysledna predpokladana uroven medianu, kvuli zaokrouhlovacim
  ! chybam by nemela byt mensi nez 10 000 pro 16bit snimky

  real, parameter :: minvalue = 0.0

  !    Step over image:
  Integer :: nstep = 1

  Logical :: simple,extend,anyf!,exist
  Integer :: i,j,n,nobs
  Integer :: istat, bpix, naxis, naxes(DIM)
  integer :: bitpix = -32
  Integer :: BlockSize,pcount,gcount

  ! image dimensions
  Integer ::  nx = -1, ny = -1 
  ! image
  Real :: skymod, skysig, medsig, medmed 
  Real :: rl, r, e, mainsig, maxvalue
  Real, allocatable :: flat(:,:), fbuf(:)
  real, dimension(:,:), pointer :: ccd
  Character(len=32768) :: sb,key,value, line, record
  Character(len=80) :: buf,dateobs,f,dateobs1,filter1
  character(len=8) :: FILTER='FILTER',DATE_OBS='DATE-OBS',EXPTIME='EXPTIME'
  type(Imagetype), pointer :: list,curr


  call InitList(list)
  curr => list
  medmed = huge(medmed)

  do
     read(*,'(a)',iostat=i,end=20) record
     if( i /= 0 ) stop 'Failed to read of an input.'
  
     call input_rcomment(record,line)

     select case(input_classify(line))

     case(input_key)

        call input_keys(line,key,value)

        if( key == 'Output filename' ) then
           read(value,*,iostat=i) nameflat
        endif

        if( key == 'Output level' ) then
           read(value,*,iostat=i) level
        endif
             
        if( key == 'Grid step' ) then
           read(value,*,iostat=i) nstep
        endif

        if( key == 'Output BITPIX' ) then
           read(value,*,iostat=i) bitpix
        endif

        if( key == 'Key FILTER' ) then
           read(value,*,iostat=i) FILTER
        endif

        if( key == 'Key DATE-OBS' ) then
           read(value,*,iostat=i) DATE_OBS
        endif

        if( key == '' .or. value == '' .or. i /= 0 ) then
           write(*,*) 'This record: ',trim(key),trim(value),' parse failed.'
        endif

     case (input_data)

        read(line,*,iostat=i) sb
        if( i /= 0 ) then
           write(*,*) 'Failed to read of filename:',trim(line)
           stop 'Bad format of the input filename.'
        end if
        write(*,*) trim(sb),':'

     istat = 0
     ! input initialisation of a current file
     Call ftopen(25,sb,0,BlockSize,istat)
     Call ftghpr(25,DIM,simple,bpix,naxis,naxes,pcount,gcount,extend,istat)
     if( istat /= 0 ) goto 666

     if( getno(curr) > 0 .and. (nx /= naxes(1) .or. ny /= naxes(2)) ) then
        write(*,*) "The image dimensions does not match. Skipping."
        call FTCLOS(25,istat)
        goto 665
     endif

     dateobs = ''
     CALL ftgkys(25,DATE_OBS,dateobs,buf,istat)
     if( istat /= 0 ) then
        write(*,*) 'DATE-OBS keyword not found.'
        istat = 0
     endif

     e = 0.0
     CALL ftgkye(25,EXPTIME,e,buf,istat)
     if( istat /= 0 ) then
        write(*,*) EXPTIME,' keyword not found.'
        istat = 0
     endif

     f = ''
     CALL ftgkys(25,FILTER,f,buf,istat)
     if( istat /= 0 ) then
        write(*,*) 'FILTER keyword not found.'
        istat = 0
     endif

     ! image - data input
     allocate(ccd(naxes(1),naxes(2)),stat=i) 
     if( i /= 0 ) then
        write(*,*) 'Not enought memory. Skipping.'
        call FTCLOS(25,istat)
        goto 665
     endif
     call FTG2DE(25,1,minvalue,naxes(1),naxes(1),naxes(2),ccd,anyf,istat)
     call FTCLOS(25,istat)
     if( istat /= 0 ) goto 666

     ! first image => initialisation
     if( getno(curr) == 0 )then
        nx = naxes(1)
        ny = naxes(2)

        Allocate(fbuf(nx*ny),stat=i)
        if( i /= 0 ) Stop 'Not enought memory.'
        fbuf = 0.0
        dateobs1 = dateobs
        filter1 = f
     endif

     ! z kazdeho nstep-ho pixelu obrazku se udela lokalni median
     call linearize(ccd,nstep,fbuf,n)
     call rmean(fbuf(1:n),skymod,skysig)
     write(*,*) '  No., Median, skysig, n:',n,skymod,skysig,n

     if( abs(skymod - level) < abs(medmed - level) )then
        medmed = skymod
        medsig = skysig
     endif

     ! add image
     Call AddItem (curr,array=ccd,skymod=skymod,skysig=skysig, &
          dateobs=dateobs,exptime=e,filter=f)

665  continue

     end select

  enddo

  ! input of top line with init values
!  read(*,*,end=20,iostat=i) sb,e,bpix
!  if( i /= 0 ) stop 'Bad format of input data.'
!  if( sb /= '' ) nameflat = sb
!  if( e > 0 ) level = e
!  if( bpix /= 0 ) bitpix = bpix

!  do
!     ! nacteni jmena souboru, ze STDIN
!     read(*,*,end=20,iostat=i) sb
!     if( i /= 0 ) stop 'Bad format of input filename.'
!     write(*,*) trim(sb),':'
!
!  enddo

20 continue
  nobs = getno(curr)
!  if( bitpix == 0 ) bitpix = bpix
  
  if(  nobs == 0 ) stop 'No input image.'

  write(*,*) ' Number of input images: ',nobs

  ! prumerna hodnota SKYSIG na urovni 'level'
  mainsig = level/medmed*medsig
  write(*,*) 'mainsig,medmed,medsig=',mainsig, medmed, medsig 

  ! prepocet vsech obrazku tak, jako by mely stejne SKYSIG
  curr => list
  do 
     curr => GetNext(curr)
     if( .not. associated(curr) ) exit

     call GetItem(curr,array=ccd,skymod=skymod,skysig=skysig)
     ! intenzity se nasobi timto cislem, aby se dostala intenzita 'level'
     rl = level / skymod
     ! prepocet 'SKISIG' pro tento obrazek na standardni
     r = mainsig / skysig
 
     write(*,*) getno(curr),' : ',r,rl,' (level, sig)'

     !        ccd = rl*ccd
     ccd = r*(ccd - skymod) + level

  enddo

  ! konecny vypocet flat-fieldu
  Allocate(flat(nx,ny),stat=i)
  if( i /= 0 ) Stop 'Not enought memory.'
  Do j = 1,ny
     Do i = 1,nx

        curr => list
        do 
           curr => GetNext(curr)
           if( .not. associated(curr) ) exit
           call GetItem(curr,array=ccd)
           n = getno(curr)
           fbuf(n) = ccd(i,j)
        enddo
           
        call rmean(fbuf(1:nobs),skymod,skysig)

        if( bitpix > 0 ) then
           maxvalue = 2**bitpix - 1
           if( skymod < minvalue ) skymod = minvalue
           if( skymod > maxvalue ) skymod = maxvalue
        end if
           
        flat(i,j) = skymod
        
     enddo
  enddo

  ! konecny vypocet medianu
  call linearize(ccd,nstep,fbuf,n)
  call rmean(fbuf(1:n),skymod,skysig)

  write(*,*) '     Final median, skysig, n=',SKYMOD,skysig,n

  ! remove a previous autoflat image
!  Inquire(file=nameflat,exist=exist)
!  if( exist ) then
!     open(26,file=nameflat)
!     close(26,status='DELETE')
!  endif

  ! cut-off range 
  if( bitpix > 0 )then
     maxvalue = 2.0**bitpix - 1.0
     where( FLAT < minvalue )
        FLAT = minvalue
     endwhere
     where( FLAT > maxvalue )
        FLAT = maxvalue
     endwhere
  end if
  
  ! Vytvoreni vystupniho obrazku
  call ftinit(26,nameflat,1,istat)
  ! zapsani informace do hlavicky
  ! dimenze a ostatni parametry se berou podle posledniho obrazku.
  ! pokud maji obrazky pretim odlisne dimenze, vysledek je nedefinovan 
  extend = .true.
  call ftphps(26,bitpix,naxis,naxes,istat)
  call ftpkys(26,'OBJECT','Flat-field',' ',istat)
  call ftpkye(26,'SKYMOD',skymod,10,'Mean level',istat)
  call ftpkye(26,'SKYSIG',skysig,10,'Standard deviation',istat)
  call ftpkys(26,FILTER,filter1,'filter of the first on input',istat)
  call ftpkys(26,DATE_OBS,dateobs1,'UTC of the first on input',istat)
  if( bitpix > 0 ) then
     call ftpkye(26,'BSCALE',1.0,10,'',istat)
     call ftpkye(26,'BZERO',2.0**(bitpix-1),10,'',istat)
  endif
  
  if( nobs > 0 ) then
     write(buf,'(a,i0,a)',iostat=i) 'The image is a result of robust average and scaling of ',nobs,' exposure(s).'
     if( i == 0 ) call ftpcom(26,buf,istat)
     call ftpcom(26,'The time of start, filter and exposure time for each image used:',istat)
     curr => list
     do 
        curr => GetNext(curr)
        if( .not. associated(curr) ) exit
        call GetItem(curr,dateobs=dateobs,exptime=e,filter=f)
        write(buf,'(a, 2x, a, 1x, f0.1)',iostat=i) trim(dateobs),trim(f),e
        if( istat /=0 ) buf = trim(dateobs)//' ..'
        call ftpcom(26,buf,istat)
     enddo
     call ftpcom(26,'All times are derived from DATE-OBS, FILTER and EXPOSURE keywords.',istat)
  endif
  call FTPCOM(26,version,istat)
  
  ! zapsani  dat
  call FTP2DE(26,1,nx,nx,ny,flat,istat)
  ! Zavreni souboru
  call FTCLOS(26,istat)

666 continue

  if( allocated(fbuf) ) deallocate(fbuf)
  call DestroyList(list)

  if( istat /= 0 )then
     write(*,*) 'The following error(s) has been occured:'
     do
        Call ftgmsg(buf)
        if( buf == ' ' ) exit
        write(*,*) trim(buf)
     enddo
     
  endif

contains

  subroutine linearize(ccd,nstep,fbuf,n)

    implicit none
    integer, intent(in) :: nstep
    integer,intent(out) :: n
    real, dimension(:,:),intent(in) :: ccd
    real, dimension(:),intent(out) :: fbuf
    integer :: i,j,nb
    
    nb = size(fbuf)
    n = 0
    do i = 1,size(ccd,1),nstep
       do j = 1,size(ccd,2),nstep
          if( n == nb ) exit
          n = n + 1
          fbuf(n) = ccd(i,j)
       enddo
    end do

  end subroutine linearize

end Program AutoFlat

